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Abstract 

We use Monte Carlo techniques and analytical methods to study the phase 
diagram of multicomponent Widom-Rowlinson models on a square lattice: 
there are M species all with the same fugacity z and a nearest neighbor hard 
core exclusion between unlike particles. Simulations show that for M between 
two and six there is a direct transition from the gas phase at z < Zd{M) to 
a demixed phase consisting mostly of one species at z > Zd(M) while for 
M > 7 there is an intermediate "crystal phase" for z lying between z c (M) 
and Zd{M). In this phase, which is driven by entropy, particles, independent 
of species, preferentially occupy one of the sublattices, i.e. spatial symmetry 
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but not particle symmetry is broken. The transition at Zd(M) appears to be 

first order for M > 5 putting it in the Potts model universality class. For 

large M the transition between the crystalline and demixed phase at Zd(M) 

can be proven to be first order with Zd{M) ~M— 2 + 1/M+..., while z c (M) is 

argued to behave as X cr /M, with A cr the value of the fugacity at which the one 

component hard square lattice gas has a transition, and to be always of the 

Ising type. Explicit calculations for the Bethe lattice with the coordination 

number q = 4 give results similar to those for the square lattice except that 

the transition at Zd(M) becomes first order at M > 2. This happens for all 

q, consistent with the model being in the Potts universality class. 
PACS numbers: 64.60.Cn, 05.50.+q, 02.70.Lq, 75.10.Hk 
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I. INTRODUCTION 



In 1970 Widom and Rowlinson (WR) introduced an ingeniously simple model for the 
study of phase transitions in continuum fluids [[I]. It consists of two species of particles, 
A and B, in which the only interaction is a hard core exclusion between particles of unlike 
species, i.e. the pair potential v a p(r) is infinite if a ^ f3, and r < Rab, and is zero otherwise. 
WR showed how this model can be transformed (by integrating over the coordinates of 
one species) into a one component model with explicit many body interactions. The A- 
B symmetry of the demixing phase transition in the original model, assumed by WR to 
occur in dimensions v > 2 when the fugacity, za = zb = z is large, then yields interesting 
information about the corresponding liquid- vapor transition in the transformed model Jl|. 

A rigorous proof of the existence of a demixing transition in this model was given by 
Ruelle 0. Ruelle used a brilliant adaptation of the Peierls argument for the Ising model on a 
lattice which exploits the A-B symmetry. Further results were obtained in ||. Ruelle's proof, 
which permits also a smaller hard core Raa = Rbb < (v / 3/2)-Rab between like particles, 
was generalized by Lebowitz and Lieb J|] to the case where vab(j) is large positive but not 
infinite. An extension of the proof to non-symmetric multicomponent models was made by 
Bricmont, Kuroda and Lebowitz using Pirogov-Sinai theory ||. We refer the reader to || 
and references there for additional results on these WR models which are, as far as we know, 
the only continuum systems with fixed decaying potentials where one has been able to prove 
rigorously the existence of phase transitions. 

The lattice version of the multicomponent WR model — hard core exclusion between 
particles of M different species on nearest neighbor sites of a simple cubic lattice in v- 
dimensions — is much easier to handle rigorously. A proof of the demixing transition and 
much more, e.g. existence of sharp interfaces between coexisting phases, in v > 3, at large 
fugacity z > Zd(M), can be obtained using standard Peierls methods, see ||, [0. A rather 
surprising result (at least on first sight) was found by Runnels and Lebowitz ||. They 
proved that when the number of components M is larger than some minimum M then the 
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transition from the gas phase at small values of z to the demixed phase at large values of 
z does not take place directly. Instead there is, at intermediate values of z, z c < z < z d , 
an ordered phase in which one of the sublattices (even or odd) is preferentially occupied, 
i.e. there is a crystalline (antiferromagnetically ordered) phase in which the average particle 
density on the even and odd sublattices, p e and p Q are unequal. The average density of 
species I = 1, . . . , M, p(I), is the same on each sublattice and equal p e (I) = M _1 p e and 
p {I) = M _1 p Q . The nature of the symmetry breaking is thus very different from that in the 
demixed phase at z > z d where p e = p D = p but there exists one species, say I', for which 
p(I') > M _1 p. The origin of the spatial symmetry breaking leading to the crystal phase 
is purely entropic. For z fixed, and M going to infinity it pays for the system entropywise 
to occupy just one sublattice without any constraint; there being no interaction between 
particles on the same sublattice each site can be occupied independently by a particle of any 
species, i.e. if we keep one of the sublattices empty then there are M independent choices at 
each site of the other sublattice. This more than compensates, at M > M , for the "loss" 
of fugacity occasioned by keeping down the density in the other sublattice. Put in another 
way, for M large enough, the typical occupancy pattern on a lattice (ignoring the label / 
of the particles) should behave like a one component lattice gas with nearest neighor hard 
core exclusion for which Dobrushin || proved the existence of a crystalline state. 

Once this is understood the natural question is, just how big does M have to be to 
see this ordered phase for M > M . It was shown in || that M < 27 6 ; a ridiculously 
large upper bound. On the other hand application of Pirogov-Sinai theory as in [|Kj where 
a similar model, in which there is a positive energy U < oo when neighboring sites are 
occupied by different species, is considered can probably be made to give M ~ 25 for our 
model corresponding to U = oo. Furthermore, a direct computation on the Bethe lattice 
with g-neighbors, gives M = [q/ (q— 2)] 2 , which would suggest M ~ 4 for the square lattice. 
The Monte Carlo (MC) simulations presented here grew out from a desire to answer this 
question and to get a general picture of the phase diagram of this system in the M — z 
plane. This we succeeded in doing, with the large M analytic results smoothly matching up 
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with the MC small M results and with the "dilute lattice model" investigated in ||10|| . The 
outline of the rest of the paper is as follows. In Section II we investigate the large M behavior 
of the model, Bethe lattice computations are given in Section III and MC simulations in 
Section IV. 



II. ASYMPTOTICS OF THE PHASE DIAGRAM 

We consider a two dimensional square lattice Z 2 . Each lattice site can be either empty 
(/ = 0) or singly occupied by a particle of type I = 1, 2, • • ■ , M . All the components have the 
same activity z and there is an infinite repulsive interaction between particles of different 
type on nearest neighbor sites. Thus a particle of type / at lattice site i can only have 
vacancies or particles of type / on nearest neighbor lattice sites. The interaction potential 
<Pi,j(hj) between a particle of type I at site i and a particle of type J at site j is: 

!oo if i and ?' are nearest neighbors and J^J, 7^0, J^O 
(1) 
otherwise 

It is clear that if we replace oo in ([!]) by some U ^ then our system is equivalent 
to a dilute Potts model. In such a system some lattice sites are empty while others are 
occupied, with a weight given by the fugacity z, by an M-component Potts variable with 
nearest neighbor interaction between like states equal to —U. The WR system (|I|) can thus 
be considered as the zero temperature limit of such a model with U > 0. We refer the reader 



to [10] for a general discussion of the phase diagram of such models. 



A. The Boundary Between Disordered and Crystal Phases 

We will argue here that the asymptotics, as M — ► oo, of the boundary z = z c (M) between 
disordered and crystal phases is given by the hyperbola Mz = X cr , where A cr is the critical 
fugacity of the one-component lattice gas with the nearest-neighbor hard-core exclusion. 

In the latter model every site of the lattice can be occupied by a particle with fugacity 
A > 0. Hard-core repulsion requires that occupied sites are not nearest neighbors. It is well 
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known |TT] that for this model there exists a critical fugacity A cr ~ 3.7962 such that for 
A < A cr the model posseses a unique limit Gibbs state (disordered phase) while for A > A cr 
it has at least two different limit Gibbs state (crystal phases). In one of the crystal phases 
the probability of the even site to be occupied by a particle is greater than that of the odd 
site and vice versa for the other crystal phase. 

Consider now a model with M different species in which we forbid particles of any type 
to be nearest neighbors. Then, for fixed A = Mz, this system is equivalent to the one- 
component lattice gas just described. Hence the multi-component WR system, where parti- 
cles of the same type can occupy neighboring lattice sites is, in a sense, the one-component 
system with the hard-core condition being slightly relaxed. That naturally leads to a con- 
jecture that for Mz < A cr our multicomponent system is in the disordered phase. 

Speaking more precisely, any configuration of the multicomponent system ([]]) can be 
uniquely decomposed on the connected components of the occupied sites. The connected 
component consisting of a single site can be interpreted as a hard-core particle with the 
fugacity Mz as we do not specify the type of the particle. Other connected components 
consisting of n > 2 sites we call clusters. All the particles in the cluster are of the same type 
and the fugacity of the cluster is Mz n , as we again do not specify the type of the particles 
in the cluster. Note that this representation naturally extends the multicomponent model 
to non integer values of M. 

It is obvious that the fugacity of clusters tends to as z — > 0, A = Mz fixed, and the 
free energy of the multicomponent model tends to the free energy of the one-component 
hard-core gas with the fugacity A, i.e. the contribution of the clusters becomes negligible in 
the limit z — > 0, with A = Mz fixed. This suggests that in the (M, z) plane there exists a 
curve M = M c (z), with M c (z) — > X cr /z as z — > 0, on which there is the second order phase 
transition between disordered and crystal phases. The typical configuration of the crystal 
phase has one sublattice, say the even one, occupied by particles except for rare excitation 
"islands", where every "island" is either a cluster or has the structure of the opposite crystal 
phase with occupied odd sublattice. 
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B. The Boundary Between Crystal and Demixed Phases 



The boundary between crystal and demixed phases admits a rigorous analysis in the 
framework of the Pirogov-Sinai theory (see |12[| , ]T3| , [H|). Related general results for the 



wider class of models including our multi-component system can be found in [plj. Here we 
present a more detailed analysis by means of a direct "elementary" approach. 

Given a configuration of the multi-component system we say that this configuration at 
a site of the lattice Z 2 is, 

(i) in the demixed phase /, if this site and at least one of its nearest neighbor sites are 
occupied by particles of type J; 

(ii) in the odd crystal phase, if this site is even and empty or is odd and occupied by a 
particle (of any type) while all four nearest neighbor sites are empty; 

(iii) in the even crystal phase, if this site is odd and empty or is even and occupied by a 
particle (of any type) while all four nearest neighbor sites are empty. 

It is not hard to check that the statistical weight of an arbitrary configuration can be 
calculated as the product over all unit bonds of the lattice Z 2 of the following statistical 
weights of the bonds: 

(i) z 1 ! 2 for the bond joining two sites in the same demixed phase; 

(ii) (Mz) 1 / 4 for the bond joining two sites in the same crystal phase (necessarily one of 
the sites is occupied and other is empty); 

(iii) 1 for the bond joining sites in the different crystal phases (necessarily both sites are 
empty); 

(iv) z 1//4 for the bond joining sites in the crystal and demixed phases (the site in the 
crystal phase is necessarily empty). 

All other bonds are forbidden, i.e. they have a zero statistical weight. Without changing 
the Gibbs distribution we multiply all statistical weights of the bonds by z~ x l 2 and obtain 
renormalized statistical weights 1, (M/z) 1 ^ 4 , z~ x l 2 and z _1//4 respectively. From that picture 
it is clear that for 1 < M < z the configurations with all sites being in the same demixed 



phase are the only periodic ground states (i.e. the configurations minimizing the specific 
energy) of the model (there are M of them). For < z < M the configurations with all sites 
being in the same crystal phase are the only periodic ground states of the model (there are 
2 of them). Finally on the line M = z all M + 2 configurations above are the ground states 
of the model and there is no other periodic ground state. 

Given a configuration we introduce the Peierls contours as the connected components 
of the unit bonds of the dual lattice Z 2 = Z 2 + (1/2, 1/2) separating sites of Z 2 not in the 
same phase. It is not hard to see that every configuration has an equivalent representation 
in terms of a collection of mutually disjoint Peierls contours. Moreover, on the line M = z 
the statistical weight of the configuration is the product of the statistical weights of the 
contours. In turn, the statistical weight of a contour is the product of the statisical weights 
of the bonds of Z 2 dual to the bonds of the contour. Here dual means rotated by 7r/2 with 
respect to the center of the bond. This representation is customary for the Pirogov-Sinai 
theory and allows complete investigation of the diagram of the periodic limit Gibbs states 
in the region M > max(z /z, 1) for a sufficiently large absolute constant Zq. 

Theorem. For M > max(zo/2, 1) and Zq large enough there exists in the (M,z) plane 
a curve M = M^(z), \Md(z) — z\ — 0(1), such that on this curve all M + 2 ground states 
generate the corresponding limit Gibbs measures and these limit Gibbs measures are the 
only periodic limit Gibbs measures. 

For max(z /z, 1) < M < M^{z) every demixed ground state generates the corresponding 
limit Gibbs state and these M limit Gibbs states are the only periodic limit Gibbs measures. 

For M > max(z /z, 1, Md(z)) every crystal ground state generates the corresponding 
limit Gibbs state and these 2 limit Gibbs states are the only periodic limit Gibbs measures. 



(See 0, 0, and cf. |T(J). 

More precise asymptotics Md(z) = z + 2 — z~ l + 0(z~ 2 ) as z — > oo can be calculated 
via an approach suggested in (jOJ. Namely, set Md(z) = z + Y^=oC n z~ n . Then one can 
successively calculate c n by equating term by term the cluster or polymer expansion series 
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written for the specific free energies f c and fd of the crystal and demixed phases respectively. 
The definition of polymers and their statistical weights can be found in [|1(| . 
It is not hard to see that 

f d = -\ogz-z- l + ^z- 2 + 0{z-*), 

where 

(i) — \ogz is the specific energy of the demixed ground state; 

(ii) z~ l is the statistical weight of the polymer consisting of a single excitation obtained 
from the demixed ground state by removing a particle from a given lattice site; 

(iii) z~ 2 /2 is the statistical weight of the polymer consisting of two copies of the excitation 
defined in (ii); 

(iv) 0(z~~ 3 ) is the contribution of the rest of polymers. 

Similarly 

f c = ~ log(Mz) - \iMzy 1 - X -M{zM~ A ) + 0(z~ 3 ) 
= -togz-|z- 1 -(|-| + l) z-* + 0(z-% 

where 

(i) — [log(Mz)]/2 is the specific energy of the crystal ground state; 

(ii) (Mz)~ l is the statistical weight of the polymer consisting of a single excitation 
obtained from the crystal ground state by removing a particle from an occupied site; the 
factor 1/2 is due to the fact that in the crystal ground state only one half of the lattice 
points are occupied; 

(iii) zM~ 4 is the statistical weight of the polymer consisting of a single excitation obtained 
from the crystal ground state by placing a particle of a given type 7 at a previously non 
occupied lattice site; such a particle requires the neighboring sites to be in the same phase 
/; any of M types of particles can produce such an excitation which is reflected in the factor 
M; finally, the factor 1/2 is due to the fact that in the crystal ground state only one half of 
the lattice points is non occupied; 



(iv) 0{z 3 ) is the contribution of the rest of polymers; 

Thus equating f d and f c one obtains c = 2 and c\ = —1. Note that every c n can be 
calculated in the same way but the amount of calculation grows very fast with n. 

III. THE BETHE LATTICE COMPUTATION 

We now present an exact calculation of the phase diagram for our multicomponent WR 
model formulated on the Bethe lattice of general coordination number q. The computation 



is based on the exact "inverse solution" for simply connected lattice structures |17]], |T8 
The approach is a more transparent alternative to the usual recursion method, it replaces 
the most stable fixed point criterion by the general principle of the global minimum of the 
free energy. 

In the Bethe lattice, every vertex % is an articulation point of multiplicity q. Let us 
denote by {zi(I), I = 0,1, ... , M} with the reference 2*(0) = 1 the set of fugacities assigned 
to the corresponding particle states at site i, and by {pi(I)} the generated particle-density 
field constrained by 

M 

5>W = 1 - (2) 

1=0 

A direct method assumes given fugacities {zi(I)} and then calculates the densities {pi(I)} 
and the specific free energy via the Gibbsian grand canonical ensemble formalism. In the 
inverse method on the contrary we calculate the specific free energy and fugacities {#»(/)} as 
functions of densities {pi{I)} considered as the basic variables. The articulation character 
of vertices in the Bethe lattice then permits the topological reduction of the equilibrium 
description. In particular: 

(i) the inverse profile equation, i.e. the dependence of the fugacity field {^(J)} on the 
specified density field {pi(I)}, takes the local form |Hj 

9-1 



pM 
Pi{i) 
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II :; '-(/). (3) 
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Here and later, the superscript < i,j > refers to the model (1) on an isolated (i.e. with empty 
boundary conditions) bond joining two neighboring sites i and j. The partition function 
of this two-site model is denoted by H <lJ> . The fugacities are denoted by zf i,0> '(/) and 
zf %,3> (I). The densities are {pi(I)} and {pj(I)}- In the notation for densities the superscript 
is omitted as these densities coincide with actual ones on the Bethe lattice which we always 
assume to be specified. 

(ii) The Helmholtz free energy of a system expressed via densities, F[p] = — logS + 
Pi{I) log 2* (I) with H being the grand partition function, can be calculated as (see [fL8| ) 

F[p) = J2 F <iJ> iPi, Pj] " (? - 1) E > ( 4 ) 

<i,j> i 

where the summation over < i,j > is over every pair of nearest neighbors. The one-site and 
the two-site Helmholtz free energies are 

M M 

F%] = -logS* + 5>C01°g*(I) = EM/)log Pi (J) (5) 

7=0 7=0 

and 

A/ A/ 

J^Ift, Pil = - logH <JJ> + E AW log^(/) + E PiW log%(/)- (6) 

7=0 7=0 

In order to mimic realistically a non-simply connected structure of the same coordination 
q, it is necessary to avoid the effect of the large number of boundary sites of the Cayley 
tree. This can be done by assuming that, in the thermodynamic limit, the local properties of 
interior vertices (expressed in terms of particle densities) are equivalent, i.e. by considering 
only translation-periodic extremal Gibbs measures on the Cayley tree. Under homogeneous 
external conditions Zi(I) = z(I) for all i, the existence of a phase transition can be simply 
detected from (|^) as a bifurcation point, resp. a "jump" point for the densities, with the 
minimum principle of the free energy, determining the right equilibrium. 

A. Crystal Phase 



The crystal regime of our model is characterized by the supposed symmetry-breaking in 
particle densities on alternating sublattices. Accordingly we set Pi(I) = p\ for sites i on the 
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first sublattice and Pj(I) = p 2 for sites j on the second sublattice (1 = 1,... ,M). In the 
corresponding two-site model we suppose that zf t,3> {Q) = Zj l ' 3> (0) = 1, zf l,:)> (I) = Z\ and 
zf t,3> (I) = z% (I = 1, 2, . . . , M) since there is no preference for any component at a given 
sublattice. With this notation the two-site pair partition function 2 <M> is 

= l + M(zi + z 2 ) + Mz x z 2 . (7) 

As functions of z\ and z 2 the one-site particle densities p± and p 2 are given by 

p 1 S <ij ' > = zi(l + z 2 ) • (8a) 



P2 S<^> = ^(1 + ^)- (8b) 



^From (^), (H) we easily get E <hJ> and zi j2 as functions of {pi,p 2 }: 



M(p 1 + p 2 )-l + p 1 -p 2 + D 1 / 2 
Zl = 2(1 -M Pl ) ' (9a) 



M(p 2 + p 1 )-l + p 2 -p 1 + D 1 / 2 

Z2 = w^;) ' (9b) 



7<i,3> 



M(M - l)(pi + p 2 ) - (M - 2) + MD 1 / 2 



2(1 -Mpi)(l -Mp 2 ) 
where the plus sign of the square root of the discriminant 



(10) 



D = [1 - (M - + p 2 )Y + 4 Pl p 2 (M - 1) (11) 

is fixed by the condition z lj2 — > for p 12 — > 0. Finally, using (|]) and @ with homogeneous 
external conditions Zi(I) = z for all i and 7 = 1, . . . , M , we find 



1 - M Pl 

. Pi 



4 , (12a) 



'-(^T«- 
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Because of translation periodicity the free energy per site, f c , can be determined from 
as follows 



fc 



-| logH<«> - q -^- log [(1 - M Pl )(l - M P2 )] . 



(13) 



The relations (IT4 ) can be considered as the equation z(pi, p 2 ) = z(p 2 , pi) (0 < pi, p 2 < 
M _1 ). In the variables s = {p\+p 2 )/2 and t = p\—p 2 it can be rewritten as z(s, t) = z(s, —t) 
and always has a trivial solution (s,0). There is no other solution if s G [0, s% ) U (s^ , M _1 ], 

1 



2M 



2M 



(14a) 



(14b) 



and 



£ = 1 



AM(q- 1) 



(15) 



(M — l)q 2 ' 

because for these s, calculated from the equation dz(s, 0)/dt = 0, the derivative dz(s,t)/dt 
has a constant sign. For s G (sf, s^) one has a nontrivial solution to z(s, t) = z(s, —t). The 
corresponding critical fugacities read 



M q - 



.! f l + E^ \ 
{ 1-E 1 / 2 J 



(q-2)/q- £ 1/2 f 
1 - E 1 / 2 



-2)/g + E 1 / 2 ] g 
1 + E 1 / 2 



(16a) 



(16b) 



The critical point exists provided that E > 0, which imposes the requirement on the number 
of components M > M with the minimum value M given by 
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(17) 



It is easy to verify, using (|13"D, that the crystal phase is thermodynamically dominant (i.e. 
it has the minimal specific free energy) with respect to the disordered one in the interval of 
activities zz 1 < z < z^ . It is not hard to see that for large M 
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M {q-2)i M 

and 



c 71 /t /„ o\n i\/r V / 



B. Demixed Phase 



In the demixed phase regime, the sites are equivalent but one of components, say 1=1, 
has greater density: {pi(l) = p(l), Pi(I) = p(2) for all 7 = 2,..., M}. Assuming zf l ' 3> (0) = 
zf' J> (0) = 1 and denoting z?' j> {l) = zf' j> (l) = z(l), zf' J> (I) = zf' j> (I) = z(2), I = 
2, . . . , M for the two-site model we have 

E <iJ> = [ X + z{ iy2 + (M _ + ^ (2)] 2 _ ( M - 1) , (20) 



p(l)5<*> = z(l)[l + z(l)]. 



(21a) 



p(2)S<"> = z(2)[l + z(2)] 



(21b) 



The system of equations (||) is closed by considering Z{(I) = z, I = 1, . . . , M and (@), which 
leads to the equations 

■l-p(l)-(M-l)p(2)^ ; - 1 



p(i; 



(22a) 



19-1 



1 - - (M — l)p (2) 
P(2) 



The free energy per site is readily obtained in the form 



z{2f . 



(22b) 



Id = — lo 



9 i _ _ ^<i,j> 



1) log[l - p(l) - (M - l)p(2)] . 



(23) 
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We solved the non linear equations (p0|)-(p2[) numerically. Given M, there is for small 
z only one solution corresponding to the disordered phase. As we increase z, two other 
solutions appear at some z = Zd(M). For bigger values of z even more solutions exist. It 
appears that among these nontrivial solutions the solution with maximal p(l) always has 
the minimal free energy. We call it demixed 1 (dl) solution. Finally we calculate the phase 
diagram for q = 4 (see Fig. 1) by comparing the free energies of disordered, crystal and 
demixed 1 phases. 

The behavior of the solutions is very similar to that of the zero-field Potts model on the 
Bethe lattice [19|J2(| (where the role of fugacity is played by the coupling constant). More 
precisely, simultaneously with dl another solution, d2, appears at the point z d (M) (see 
Fig. 2a). At Zd, the free energy of dl and d2 phases is greater than that of the disordered 
phase and Pdi(l) — PdziX) > PdisiX)- Increasing z further, these two phases split and 
Pdi(l) > Az2(l)- At z = Zd (or z = z c ), the free energy of the dl phase coincides with the one 
of the disordered (or crystal) phase, and the system exhibits a first-order phase transition, 
accompanied by a jump in densities PdisiX) — > Pdi(l)- To complete the description, we 
mention that at a larger value of z, z — z^, given by 



we have PdziX) — Pd2(2) = Pdis and the free energies of the disordered and d2 phases coincide. 
For z > Zd, one observes that Pd2(l) < Pd2(2), i.e. one particle components is paradoxically 
suppressed by the others for d2 solution which has in this region a lower free energy than 
that of disordered phase but greater than that of dl phase. The only exception from the 
above scenario is represented by the M = 2 component WR model (Fig. 2b). In that case, 
the dl and d2 phases are in fact the equivalent realizations of the particle 1 <-> 2 exchange 
symmetry of the same demixed phase, and the corresponding demixing phase transition is 
of second order. 

Let M c [y) denotes the "critical" number of components of the WR model in v- 
dimensions, such that the phase transition from the disordered to the demixed phase is 




(24) 
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second order for M < M c and first order for M > M c . For the Bethe lattice, we have 
M c = 2 independently of the coordination number. This value of M c is exactly the same 
as the value of its counterpart defined for the ordinary zero-field M-state Potts model on 
the Bethe lattice. The equality M c = 2 for the Potts model is supposed to hold for regular 
lattices in dimensions v > 4, where the mean-field treatment provides an adequate descrip- 
tion of the critical behavior. We therefore suggest that our M-component WR model is a 
dilute version of the M-state Potts model, preserving the Zm symmetry among the particle 
states, which falls into the same universality class. This conjecture is supported by the MC 
estimate M c = 4 for the WR model on the v — 2 square lattice in accordance with the 



behavior of the Potts model on the square lattice |2T|. This is discussed in the next part of 
the paper. 



IV. MONTE CARLO SIMULATION 

In this section we present results of the Monte Carlo study of the M- component WR 
model on a square lattice of size S 2 = 100 x 100 with periodic boundary conditions. On 
an initially empty lattice we deposit particles chosen at random from the M components at 
fugacity z respecting the exclusion given by (|l|). We sequentially update the lattice using 
a checkerboard algorithm resulting in a good vectorization. An update of a lattice site 
occupied by a particle of type I (I = indicating an empty site) is done as follows: 
We randomly choose a new trial particle of type 7j r , where It r can have any integer value 
between and M with equal probability. It r = refers to a removal attempt of a particle 
1^0 from the lattice site, which is successful, if a number X randomly chosen with equal 
probability between and 1 is smaller than the inverse fugacity 1/z. In this case I gets 
the value 0, otherwise it remains unchanged. It r ^ refers to a deposition attempt of a 
particle of type Jj r . If / = then it is successful if each of the four nearest neighbor sites is 
either empty or occupied by a particle of the same type (It r ) and X < z. In this case / gets 
the value I tr , otherwise it remains unchanged. A direct replacement attempt of a particle 
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I 7^ surrounded by four empty nearest neighbor sites is always successful. Typically in 
a simulation run after an equilibration of 5 x 10 5 Monte Carlo steps (MCS) we update the 
lattice 5 x 10 5 times (in the cases M = 6 and z = 2.5, 3., 3.5 and 4 up to 5 x 10 6 times), 
the configuration of every tenth step is taken for the evaluation of the averages. A typical 
run with 5 x 10 5 MCS took about 3 CPU hours on a CRAY-YMP. 

Let m(«i,i 2 ) denote the occupancy of a site, m(ii,z 2 ) = if the site (11,12) is empty and 
m(ii,i2) = 1 otherwise. As observables we took histograms Pl{4> c ) of the order parameter 
(f) c for the crystal structure and Pii'Pd) of the order parameter (pd for the demixed phase in 
subsystems of size L x L, 

^c=^ [MMa) - 1] (-l) 414 * (25) 

and 

d = -iMa X/ Ar z (I)-p/M (26) 

where Nl(I) denotes the number of particles of type J in a subsystem of size L x L and p 
is the average overall density. 

In Fig. 3 we show typical configurations for M = 9 at three different values of z. In 
Fig. 3(a), z = 0.1, in this case the system is in the gas (or disordered phase). In Fig 1(b), 
z — 5, and the system is in the crystal phase where one of the sublattices has a higher 
density. In Fig. 3(c), z = 8.5 and the lattice is predominantly occupied by particles of one 
type (demixed regime). 

We first discuss the phase transition from the gas to the crystal phase. For a given M 
the transition activity z c is found by finite size scaling techniques In particular the 

k — th moments of the order parameter distribution Pl(4>c], 

< (f> k c > L := J <p k c P L (<p c )d<p c (27) 

can be evaluated in subsystems of size L x L, and from them the fourth order cumulant EB 



Ul , 
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In a one phase region far away from a critical point the subsystem size typically can be 
chosen larger than the correlation length £, L » £ and the order parameter distribution is 
to a good approximation a Gaussian centered around 0, resulting in Ul — > for L — > oo. In 
the two phase coexistence region far away from a critical point we can again assume L » £ 
and the order parameter distribution is bimodal resulting in Ul — ► 2/3 for L —>■ oo. Near 
the critical point however we have L << £, and using scaling arguments |23] the cumulant 



is a function of resulting for £ — > oo in the same value of U* for all different L. This 
method efficiently allows the localization of critical points by analyzing the cumulants for 
different values of z on different length scales L. For low values of z we are in the disordered 
one phase region, here Uy > Ul for U < L. For large enough M we obtain a crystal phase 
with Ul> < Ul for L' < L. Near a critical point we should expect Uy ~ Ul for U ^ L. 

In Figs. 4-6 we present the results for Ul as a function of z. For M = 9 we observe in 
Fig. 4 a cumulant intersection near z c = 0.85 ±0.05 indicating the phase transition from the 
disordered to the crystal phase. For M = 8 we have an intersection point at z c = 1.1 ± 0.05 
(see Fig. 5) and for M = 7 at z c = 1.6±0.1 (see Fig. 6). For M = 20, 15, 14 and 10 we obtain 
z c = 0.24, 0.35, 0.38, 0.68 ± 0.01 respectively. The transition points are presented in Fig. 7 
together with the asymptotic expression |24| Mz = 3.7962, (section II) and the results from 
the Bethe lattice computation, (section III). We find good agreement of our MC results with 
those of the asymptotic form for large values of M and for large z with the results of the 
Bethe lattice computation as well. 

The case M = 6 is analyzed in Fig. 8. In all cases studied the order parameter cumulant is 
decreasing with increasing system size, see Fig. 8(a), indicating the presence of the disordered 
phase. In Fig. 8(b) we show the cumulant for a given system size versus z, no cumulant 
intersection points were found in these cases, but Ul> > Ul for V < L even for z- values 
as large as z = 4.4 indicating that the system with M = 6 components is in the disordered 
phase. Our conclusion of the data analysis is that, based upon the present statistical effort, 



no sign of a phase transition from the gas phase to the crystal phase was found for M < 6, 
indicating the nonexistence of this transition for M < 6. Thus the minimum number 
of components required for the existence of the crystal phase is M = 7; with possibly a 
noninteger value M between 6 and 7. 

The transitions from the disordered to the crystal phase were analyzed further by finite 
size scaling techniques. We plotted the scaling functions of the order parameter and the 
order parameter susceptibility, 4> c = L 13 ^ < |0 C | > L and x~ c = L 2-7 ^ [< 4> 2 C > L — < \<f> c \ > 2 ] 
versus the scaling variable t = \z — z c \L l ^ v , where (3, 7 and v are the critical exponents of the 
order parameter, susceptibility and correlation length respectively. On a double logarithmic 
plot we obtained data collapses and the 2D-Ising asymptotic large-t behavior for all cases 
studied by utilizing the critical exponents of the 2D-Ising universality class (/3 = 1/8, 7 = 7/4 
and v — 1); examples are shown in Figs. 9, 10, 11. These data, in conjunction with the 
cumulant intersection value being independent of M and close to the accepted value for the 
2D Ising class indicate that independent of the number of components M the transition 
from the disordered to the crystal phase belongs to the 2D Ising universality class. 

We now discuss the phase transitions to the demixed phase. The phase transition was 
analyzed by a study of the order parameter distribution P L ((p d ), the density p(z) and free 
energy integration. In Fig. 12 we show the density of the system versus z for different values 
of M. The density is an increasing function of z and approaches for large z the asymptotic 
form p — z/(l + z) which describes the behavior of the system with M — 1. In the demixed 
phase one particle type is dominant and so the system properties are close to that of a one 
component system. In general the density in the demixed phase exceeds the value 1/2, for 
M > 6 we observe a direct first order transition from the crystal to the demixed phase, 
with a finite jump in the density at the transition fugacity z&. For large M the density 
should jump from p m 1/2 to p ~ z d /{ \ + z d ). In the simulations we find a hysteresis region 
around z d going approximately between these two values when increasing and decreasing 
the fugacity. In cases of a small hysteresis region with extent of less than \z — z d \ < 0.1 the 
middle z-value of this region was taken as the transition value z d . In cases of a pronounced 
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hysteresis region we located the transition fugacity Zd by two independent methods, (i) We 
integrated the free energy from z — to and from z = oo to Zd, where the limiting free 
energies are known. For the low fugacity free energy we obtain: 



Pl (z) = - / dz 1 p(z 1 )/z 1 (29) 
Jo 

For the high fugacity free energy we obtain: 

p h (z) = - f dzip(zi)/zi - log(l + zi) (30) 

Jzi 

where z\ is chosen that large (typically z\ = 50) that p(zi) ~ + z{). The intersection 
value Zd, where pi(zd) = Ph(zd) was one estimate for the value of the fugacity where the first 
order transition to the demixed phase takes place, see Fig. 13. (ii) Another independent 
estimate for z& was obtained by the procedure of searching for relative stability of one of 
the two phases during a simulation starting from configurations with both phases present in 
parallel slices extending over the length of the simulation box, see Fig. 13. Both independent 
methods gave, within 5% uncertainty, the same numerical values for Zd- The resulting 
phase transition values of Zd are shown in Fig. 7. We note that with increasing number of 
components, the transition fugacities approach the exact asymptotic line M = z + 2 — 1/z, 
see section II. 

For the cases M > 5 we find a jump in the density at Zd so that we classify these 
transitions as first order transitions. For M < 4 we do not observe a jump in the density at 
Zd- In these cases the transition was located by the cumulant intersection method, see above, 
where the order parameter <pd is the order parameter of the demixed phase and <p c has to be 



replaced by <pd in Eqs. fl2~7|) and (p8|). A finite size scaling analysis shows that the data for 
the order parameter and the susceptibility of the phase transitions for M < 4 are consistent 
with the 2D M-state Potts universality class, see Figs. 14, 15. For the susceptibility and 
M = 4 we obtain a much better scaling behavior compared to an analysis assuming the 
2D-Ising exponents to be valid, however we note that the ratios for (3/v and 7/1/ of the 



Potts classes |21j agree with the 2D Ising class values within a few per cents, so that a 
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clear distinction between these classes based on our numerical data is difficult. Recently the 
case M = 2 was studied with Monte Carlo methods where evidence for the 2D-Ising 
universality class was found, in agreement with our findings. 
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FIGURES 

FIG. 1. Phase diagram in the M-z plane for a square lattice (MC) and for a Bethe lattice 
with coordination number q = 4 . Lines: Exact results for the Bethe lattice for the transition lines 
from the gas phase to the crystal phase (dashed line), from the gas to the demixed phase (full line) 
and from the crystal to the demixed phase (dotted line). Symbols for MC: Transition points from 
the gas phase to the crystal phase (circles), from the gas to the demixed phase (triangles) and from 
the crystal to the demixed phase (squares). 

FIG. 2. A schematic plot of the interplay among the disordered (dis), demixed 1 (dl) and 
demixed 2 (d2) Bethe solutions in the [z, p(l)] plane for a) M > 2, b) M = 2; the solution with 
the lowest, intermediate and highest free energy is depicted by the solid, dashed and dotted line, 
respectively. 

FIG. 3. Typical configurations for M = 9 and (a) z = 0.1, (b) z = 5 and (c) z = 8.5. 

FIG. 4. Cumulant Ul versus z for M = 9. The different symbols refer to different subsystem 
sizes as indicated in the figure. Lines are for visual help. 

FIG. 5. Cumulant Ul versus z for M = 8. The different symbols refer to different subsystem 
sizes as indicated in the figure. Lines are for visual help. 

FIG. 6. Cumulant Ul versus z for M = 7. The different symbols refer to different subsystem 
sizes as indicated in the figure. Lines are for visual help. 

FIG. 7. Phase diagram in the M-z plane for a square lattice (MC). Full lines: Asymptotic 
lines for the phase transitions in the high fugacity region, M = z + 2 — l/z, and for the transition in 
the low fugacity region, Mz = 3.7962. Dashed lines: Results from the computation on the Bethe 
lattice with coordination number q = 4. Symbols for MC: Transition points from the gas phase to 
the crystal phase (circles), from the gas to the demixed phase (triangles) and from the crystal to 
the demixed phase (squares). 
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FIG. 8. (a) Cumulant Ul versus inverse system size for M = 6 and various values of z. (b) 
Cumulant Ul versus z for M = 6. The different symbols refer to different subsystem sizes as 
indicated in the figure. Lines are for visual help. 

FIG. 9. Scaling functions of the order parameter (a) and the order parameter susceptibility 
(b) for M = 20 utilizing the 2D-Ising critical exponents. Lines indicate the asymptotic power law- 
behavior with the 2D-Ising critical exponents (t = \z — z c \L 1 l v ). 

FIG. 10. Scaling functions of the order parameter (a) and the order parameter susceptibility 
(b) for M = 15 utilizing the 2D-Ising critical exponents. Lines indicate the asymptotic power law 
behavior with the 2D-Ising critical exponents (t = \z — z c \L 1 l v ). 

FIG. 11. Scaling functions of the order parameter (a) and the order parameter susceptibility 
(b) for M = 10 utilizing the 2D-Ising critical exponents. Lines indicate the asymptotic power law 
behavior with the 2D-Ising critical exponents (t = \z — z c \L 1 l v ). 

FIG. 12. Density of the system versus z for different number of components M, symbols refer 
to Monte Carlo results, the connecting lines are for visual help. The full (open) symbols refer 
to the MC results obtained by starting with configurations previosly obtained for lower (higher) 
fugacities. The full line equals p = z/(l + z). 

FIG. 13. Free energy pi(z) and Ph(z) from thermodynamic integration for (a) M = 15 and 
(b) M = 10. The arrow at the z-axes indicates the transition value found by the phase stability 
study, see text. 

FIG. 14. Scaling function of the order parameter for M = 2 (a), M = 3(b) and M = 4 (c) 
using the 2D M-state Potts critical exponents. Lines indicate the asymptotic power law behavior 
with the 2D-M-state Potts critical exponents (t = \z — z c \L l / u ). 
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FIG. 15. Scaling function of the order parameter susceptibility for M = 2 (a), M = 3(b) and 
M = 4 (c) using the 2D M-state Potts critical exponents. Lines indicate the asymptotic power 
law behavior with the 2D-M-state Potts critical exponents (t = \z — z c \L l / u ). 
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